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Tine present invention generally r4S!latee to apparatus and 
metiiods for acquiring seismic signal and filtering such 
5 data. More particularly, it is directed to a method of 
acquiring and processing single sensor recordings, 
specifically to derive filters to attenuate the noiise in the 
thus acquired seismic data* 

10 BACKGROUKD OF THE IKfVENTIOKr 

Seismic data is collected to analyze the subsurface of the 
Earth, and is particularly collected in connection with 
hydrocarbon exploration and production activities. Seismic 
data for analysing subsurface structures may be collected on 

15 land or over "water • In order to obtain seismic data, an 
acoustic source is used which typically consists of 
explosives or a seisinic vi'brator on land or an impulse of - 
compressed air at sea. The eeismic signals reflected Tcty the i 
various geologic layers beneath the surface of the Harth are 

20 3cnov?n a» traces and are sensed by a large ninnber, typically 
hundreds or thousands , of sensors such as geophones on land 
and hydrophones at sea. 'Hie reflected signals are recorded 
and the results are analysed to derive an indication of the 
geology in the subsurface. Such indications may then be used 

25 to assess the likelihood and location of potential 
hydrocarbon deposits. 

seismic surveys are generally conducted using one or more 
receiver lines having a plurality of receiver station 
30 locations spaced evenly along their lengths. In a two 

dimensional (2i>} survey, a single receiver line is used and 
the acoustic source is typically positioned at various 



1 



09-JfiN-2004 105 46 FROM SCR IP LAW DEPT 
57.0589 GB NP 



points in-line with the receiver line. In a three 
dimensional survey, a plurality of parallel receiver lines 
are typically used and the acoiastic source is generally 
positioned -at various points offjset from the receiver lines. 
5 While. a 2D seisitdc survey can only create a cross- sectional 
representation of the sul)Siirface, a 3D seismic survey can be 
used to develop a three dimensional representation of the 
Bxibsurf ace. 

10 Seismic data are subject to a wide variety of noise related 
problems that can and do limit its usefulness ^ Broadly 
speaking, noise found in seismic traces is either incoherent 
or coherent. Incoherent antbient noise ^ or uncorrelated 
■^white" noise, is ubiquitous and is generally greatly 
15 attenuated through the simple .expedient of stacking, 

although extremely large individual data values ("spikes'^) 
and "bad" traces often need special attention. Coherent, or 
correlated, noi^e on the other hand cannot usually be so 
readily eliminated* Some coimnon examples of coherent noise 
20 include multiple reflections^ ground roll, air waves, guided 
waves, sideswipe, cable noise and 60 hertz power line noise,' 
Among the many known approaches to attenuating noise, there 
are space-time based or transform based methods. Space-time 
based methods operate on the time series of the acq[uired 
25 data* Transform based methods operate on data transformed , 
from the space-time domain into another domain using a 
suitable transformation. The most popular of the 2~D 
transform methods is the 2-D Fourier transform (or "f-k" 
transform) . Seismic data containing noise are transformed to 
30 * the alternative domain where noise events are more coxiipactly 
represented. If the noise energy can be located and isolated 
in the transform domairi, it is removed from the transformed 

2 
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data by filtering or muting/ i-^. by setting the values in 
the region corresponding to the noise energy equal to zero 
or some other minimal value. Finally, the transformed data, 
without the noise energy, are then inverse transformed to 
5 return them to the time and offset (i.e., untransEormed or 
"x-t") domain. 

in conventional seismic data acquisition systems data are 
inherently filtered through use of -hard-wired" 
10 {electrically connected) groups of s«isors. A group or 
receiver array delivers a single output trace (the 
normalized sum or arithmetic average of the output of all 
individual sensors of the group) at the particular receiver 
station location about which the sensors are placed. The 
15 single trace is the normalized sum or arithioetic average of 
the output of all individual sensors making up the group. 
Without rurthef processing, such a two-dimensional group has 
a spectral response that can be approximated by a f re<auenoy- 
independent 2D sine function in the wavenumber or kx-ky 
2 0 domain . 

More recently, however, seismic surveys have been performed 
using single or point receiver arrays. Such surveys offer 
the potential of recording the output of individual sensors 
25 or receivers and the inherent filtering effect of the hard- 
wired group can be replaced by filters that are better 
adapted to the nature of seismic noise and preserve more of 
the seismic reflection signals. 

30 It is therefore an object of the present invention to 

provd.de methods for processing seismic data, particularly 
methods for designing and applying filters for such data. 
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StMKLaRY OF THE 
invention includes a metlaod of de terming a digital 
filter for seismic signals coirtprising the steps of; defining 
5 contraints representing a filter for preserving signals 

representing reflection and/cs^r refractions from sub-surface 
structure and suppressing noise signals in Beismic signals; 
and using an iterative process with eadi iteration 
comprising: transforming a filter obtained from a previous 
10 iteration into a transform domain; applying in said 

transform domain first constraints; inverse transforming the 
filter with the applied constraints into a sample domain; 
and applying in said sample domain 'second constraints to 
obtain an iterated filter . 

15 

The term- ssiciple domain and transform domain are arbitrary in 
the seacise that both domain are representations of the seimic 
signals using a different set of coordinates. The transform 
' describes the change from one domain to the other and the 
20 inverse transform describes the reverse direction- 

Preferably the transform is between the wavenumber or 
frequency-wavenumber domain and and the spatial or temporal- 
spatial, domain* This transform can be performed by thewell- 
known Fourier transform* 

25 

In a preferred embodiment the filter is constrained to' a 
predefined tolerance in one domain and to a predefined, 
response outside a finite region in the other. The 
predefined response is in most case very small or zero to 
30 suppress the signal outside the finite region. 



4 
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a?he metliod- is preferable applied using a non-rectaiigular 
transform which, can be described as staggered or hexagonal, 
even though the hexagonal case is only a special case of- the 
general staggered transform used- 

5 

Application of this staggered transform yields superior 
results when applied to data ac<^ired on a staggered grid. 
Xt is therefore another aspect of the ibnvention to use 
groups of receivers or single sensor seismic receivers 

10 distributed to obtain seismic measurements on a staggered or 
hexagonal grid- A staggered grid can be characterised by 
having non-diagonal elements in the matrix that translates 
the grid points onto themselves ♦ Xn a more practical' 
definition a staggered grid can be seesn as an array of 

15 points in which every second rows of of -points are shifted 

with respect to the remaining rows or as two overlapping and,^ 
shifted rectangular grids. 

The filter of the present inventions are preferably is a 
20 zero-phase finite in^ulse response {FIR) filter, ^hey have 
at least at. least two dimensions^ but are preferably full 
three-dimensional {3D) filters. 

These and other aspects of the invention will be apparent 
25 from the following detailed description of non- limitative 
examples and drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 
FIGS. lA-'D show a desired ideal wavenxmiber (kx-ky) response 
30 of a 2D x-y filter in three differ^t cross-sections; 



5 
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PIG. 2 illustrates steps to design a filter in accordance 
with, an exanitple of tbe invention; 

Fig. 3 shows the waveniiniber (kx-ky) response of a 2D ac-y 
filter designed in accordance to the method of FIG. 2 using 
a rectangular discrete fast Fourier transform; 

Fig. 4 shows the wavenumber (kx-ky) response of a 2D x-y 
filter designed in accordance to the method of FIG. 2 using 
a hexagonal discrete fast Fourier transform? and 



FIGS, 5A,B coxtipare single receiver data sets obtained by 
processing in accordance with an exaitqple of the present 
invention (FIG. 5B) with a data set obtained with . 
15 conventional hard-wired groups and conventional noise 
filtering (FIG. 5A) after differential moveout (DMD) 
stacking. 

DBTAIIxED DESCRIPTION 
20 A filter in accordance with the present invention is 

designed to remove at least part of the noise from the data 
acquired through a seismic survey. The known noise 
characteristics and seismic signal spectrum depend on the 
various parameter such as source and receiver locations, the 
25 acoustic properties of near surface layers and many other 

parameters within or outside the control of an operator. Any 
such pre-established knowledge of the noise characteristics 
can be used to establish a set of filter parameters that 
when combined provide a description of the filter suitable 
30 to be used as an input to a data processing machine or 
coicputer . 
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a?he exatople of FIG. 1 sWs the ideal wavenumber ^esponne of 
a 2-r filter. Two schematic perpendicular cross-sections are 
BUovm in FIG. lA illustrating parameter for a functional 
description of the design in the vravenumher or {f - kx - ky) 
domain. PIG. IB shows the £ - toe cros^-section for ky - 0, 
PIG. IC shows the middle panel ehows the f - ky cross- 
section for kx = 0, and FIG. ID shows the kx - 3cy cross- 
section at f = 10 Hz. 

m the 2-D smrvey of this example, most of the reflection 
energy is expected to be incident iri, or near the in-line 
plane, vrhereas most of the energy incident in the cross-line 
direction should be' scattered ground-roll and other types of 
noise. The passband of the initial 2-D x - y filter is 
therefore designed with elliptical contours, with greater 
attenuation in the cross-line direction. The larger passband; 
in the in-line direction is designed to preserve signal ; 
cons^onents with maxlittiam wav-enuiriber and frequency content. 

The filter of FIG. 1 is characterized by three frequency 
bands separated by the frequencies fbreakl and f break2 . 
Within each frequency band there is defined a pass zone (the 
central zone) , a transition zone and an outer redact or stop 
zone. The boundaries of the zone are marked by piecewise 
linear function defined through' the slowness paraineters si 
to s6. . 

Further pa^rameters characterizing the filter of the example 
is the aspect ratio of the elliptical pass zone and the 
desired attenuation in the stop band and "the width of the 
transition zone. 
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Other parameter define the dimensions ( length r. width) of the 
filter or the filter support for a finite impulse response 
filter (FIR) as used iii this example - 

5 Filters can be characterised by parameters others than , those 
descritoed above. But independent of the specific form of 
representation^, it is always possible to define a 
representation of a desired filter in machine readable form. 

10 Each of the parameters of such an ideal filter can be 
regarded as constraints in the time domain and in the 
frequency domain and the probleocn of designing a real filter 
is an approximation optimized with those constraints* 

15 It is known that if the constraints define convex sets in 
the set of square sutomable sequences , then applying 
alternating orthogonal projections onto these would converge 
to the optimal solution- This is called the method of 
Alternating Projections onto Convex Sets (APOCS) . a?he APOCS 

20 method has been used in various technical fields outside 
seismic processing. For details of the APOCS method 
reference is made to Cetin et al . {IEEE Signal Procasing 
Magazine, vol. 14 (1997), pp-60-64 or.Proc. Int, Syn^- 
Circuits and Systems, 19S7, pp. 1044-1047) who adapted the 

25 APOCS method to the design of zero-phase 2~D finite im^julse 
response (FIR) filters- 

In the APOCS method the filter design problCTa is formulated 
30 to alternately satisfy the frequency domain constraints on 

the magnitude response bounds and sample (space and/or time) 
domain constraints on the impulse response support. The 
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20 



algorithm is iterative and each iteration requires two 2-D 
FFT corttputations • Cetin a J . (see above ) have shovm that if 
convergence is achieved, this approach results in filters 
with approximately ecjairipple response. 

It can be proven that the APOCS method is globally- 
convergent under, certain conditions. If the constraints in 
the sainple and spectral domains define convex, sets in the 
set of square summable seqiaenceer then the imposition of the 
constraints in the sanip;Le and spectral domains are the 
orthogonal projections onto these eetfe, and if the sets 
intersect, the iterates converge to a memiber in the 
intersection set* Furthermore/ if there is just one sequence 
satisfying both conditions^ then this sequence is the 
equiripple solution of the filter design problem. 
Consequently, the iterations converge to the equiripple ,^ 
filter. i 

Generalizing to the design of a filter, the filter i& to: 

be zero-phase, and its spectral response in the f-kx-fcsf 
domainr HCf^kx^ky) to be within given tolerance ranges: 



where HiaCfrkx.ky) is the ideal filter response^ and 
^di^r'kx,lky) is the tolerance, or the desired laaxiimiin ripple 
level, which sriay take different values in -different 



[1] 



Hid (f * kx ' ky) Ea (f - ' ky) S Hid (f .kx - ky) 

^ Hid(f.kj5iky) + Ed(fr3Cx/ky)r 



passbanda and etopbanda. ^e ideal response is defined as 



30 



123 
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and the tolerances as 



[3] 



E<3 



where Fp and F3 are tlxe pass and stop bands, and 5p and 5s are 
tiae corresponding tolerances, respectively. 

In the t-x-y domain^ the filter has a finite- extent Biipport^ 
10 I, which is a Si>inmetric region around the origin in order to 
have a izero-phase response. The time^space domain constraint 
requires that the filter coefficients are zero outside the 
resrion Z» 

15 The iterative method is initialized with an arbitrary 

finite-extent, real 3-D sequence that is syometric around 
the origin in order to have a zero-phase response. The 
time/space domain constraint requires that the filter 
coefficients must be equ.al to zero outside the region I- As 

20 an initial estimate, the inverse Fourier Transform (F"^) of 
the ideal frequency response is used; 



C43 



25 with 



[5] 



ho (t^ Xiry) = 



hid (t^rXr y)# if (t^3C,y) « 
0# otherwise- 



30 



FIG. 2 illustrated the iterative process beginning with the 
selection of an initial time domain response hid- 
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At each iteration of the algorithm^ the time/ space and 
frequency /waven-ujober domain constraints ars successively 
instposed onto the aurrent iterate- The n-th iteration 
consists of the following steps J 

- Compute the 3D Fourier Transform of the n-th iterate 
hn(trX,y)*by a 3D FFT algorithm to obtain Hn(f,kx,ky) - 

- Iir^ose the frequency /wavenuihber domain constraints as: 

Hid(f/3«bcrky) + Ed(f.kx'ky) 

if Ha(f.3Cxrky) > Hid (f / ' l^y ) + Ed (f , , ky) 

Hid fcrk^ ' ky) - Ed (f .kjc r ky\ 

if Hn (f / kjc . ky) < Hid {f . kjc * ky) + fia (f / k^c ' 
{-g] Hn (f ^kx '-ky)' Qt-herwise. 



- Compute the 3-D inverse Fourier Transform of Gii(f,kx,ky) 
to obtain gn(trX/y)*. 

- equal gii(t,x,y) to zero outside the region I to obtain 
hn.+i(t,X/y) as 



15 



[71 



hn+l (tr 3<> y) = 



'Sn »;f y)r if (Ux^y) e Ij 
Of otherwise/ 



- If the mean-squared error between the iterates ha(t,5C,y) 
and hn+i(t,x,y) is less than a predefined threshold, then 
exit. 



20 
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These steps are performed using a standard, numerical 
programme such as MatlabCTM) and the routines embedded 
therein. 

5 Hie method oiutlined above can be directly applied to geismie 
data acquired on a rectangular grid by using a 2D 
rectangular fast Fourier transformation (FFT) that samples 
the data at rectangularly distributed points. 

10 However, the advent of single sensor recording makes it more 
feasible to acquire seismic data on. a staggered layout of 
receivers. It is known that circularly bandlimited signals 
require 13 percent less samples to achieve the same spectral 
resolution when sampled on a staggered grid than being 

15 sampled on a rectangular grid, Since seismic signals and 
noise can be assumed to be circularly bandlimited in the 
waveriuiriber space ^ hexagonal sampling is more efficient than 
rectangular sampling potentially reducing the number of 
receivers in a suxrvey by 13 percent. 

20 

Xt is therefore advantageous to adapt "the present method to 
* - accommodate surveys performed using a hexagonal distribution 
of receivers or groups of receivers - 

25 It is feasible to apply the above method to data acquired on 
a staggered or hexagonal grid by using the following steps; 
-spatial oversampling : 
-performance of a - rectangular FFT; 

- imposing spectral -constraints in the wavenumtoer domain; 
30 - using a rectangular inverse Fourier transform to revert to 
the time domain; and* 
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- iinposingr zero^ on points tiaat do not coafrespond to 
receiver locations » 

FIG, 3 sliows the resulting reisjponse of a 2D filter witlx dB 
5 contour levels for a hexagonal receiver grid using a 

rectangular FPT as outlined in the previous paragraph- As in 
FIG. 1-, the top panel* shows the f - kx cross- section for ky 
= 0, the middle panel shows the f - cross-section for 

- Or ^nd the bottom panel shows the kx - cross-section at 
10 f = 10 Hz. Comparison with PTG. 1 shows that signals from 

outside the desired region (sidelobes) leak into the 
filtered data. 

This is a conse«2uence of a mismatch between the true 
15 receiver locations which for this exatnple were placed on a 
hexagonal grid and the rectangular DFT used for the 
iterative process ^ "liie method introduces holes into the 
sample doiciaizi, which is no longer convex. Convergence is 
then no longer guaranteed and the above filter design 
20 method is not expected to give optiiaum results - 

It was found that the performance of the filter design 
method can. be significantly improved by using a Fourier 
^ transform with hexagonal sampling in the space (x-y) domain 
25 when the data is acquired on a hexagonal grid. 

A discrete Fourier transform (DFT) that relates a hexagonal 
sampled signal to hexagonal sample of its Fourier transform 
is known as such and is given for example by 

30 
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Xt ki, k2) = 2j ^ ^1' ^2 ^ ®^ 




iil=0 n2==0 



5 For this and further detallB of the liexag-onal DFT reference 
can be made for exa^n^le to R* M, Mersereau and D, B Dudgeon 
in ^'Multidimensional Digital Signal Processing'', Prentice- 
Hall 1984, 97-100. 

10 FIG, 4 shows tlie response resulting from the use of a 

hexagonal DTF of the 2D filter with dB contour levels. As in 
FIG* If the top panel shows the f - kx cross-section for ky 
- 0, the middle panel shows the f - ky cross-section for kx 
= 0, and the bottom panel shows the kx ~ ky cross-section at 

15 f = 10 Hz. Coittparison with FIG. 1 shows a very good match 
between the desired and the resulting responses, and the 
sidelobe levels of the stop band are less than 35 dB 
compared to the higher sidelobes apparent in FIG- 3. 

20 It was foimd that the filter design results in improved data 
processing^ particularly when used on seismic data acquired 
through single sensor surveys. To illustrate the irnprovement 
the. same data set was acquired using (a) conventional groups 
of . receivers, (b) single sensors and (c) a very low density' 
25 of single sensors. After filtering the -data was further 

processed using a stacking procedure with dif f ereaatial move- 
out (DM0) correction. 

FIG. 5A shows the DMO stack obtained using conventional 
3 0 filtering methods on signals from analogue groups of 72 

receivers par 30 metres, while FIG, 5B depicts the DMO stack 
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obtained after processing by an adaptive pre-filtex and an 
2D filter desigioted" tJirough the iterative process described 
above (FIG, 2), using 36 receivers per 30 metres. Coirpared 
to the analogue result, FIG. 5B demonstrates considerable 
5 increase in S/N, especially for higher CMP's, 
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CLAIMS 



1- 



Metliod of determine a digital filter* for seismic 



signals coiaprising the steps of: 
5 defining- odrxtrainta representingr a filter for preserving- 

signals repres^ting reflection and/ or refractions from sub- 
surface structure and suppressing noise signals in seismic 
signals; and 

using an iterative process with, each iteration further 
10 comprising the steps ofi 

- transforming a filter obtained from a previous iteration 
into a transform domain; 

- applying in said transform domain first constraints ; 

- inverse transforming the filter with the applied 
15 constraints into a sample domain; and 

- applying in said saiuple domain second constraints to 
obtain aix iterated filter * 

2. The method of claim 1 wherein each step of the 

20 iterative process includes the transform of the filter 

(coefficients) into the wavenimtoer or frequency-wavenumber 
domain and the inverse transform back into the -spatial or 
temporal -spatial domain. 

25 3. ^e method of claim 2 wherein in each step of the 

iterative process the filter is constrained to a predefined 
tolerance in the wavexiuinber or f reqiuency-wavenuxnber domain. 

4- The method of claim 2 wherein in each step of the 

30 iterative process the filter is constrained to a predefined 
response outside a finite region in the spatial or temporal - 
spatial domain, . 
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5. The method of claim 2 wherein in eacli step of tlie 
iterative process the filter is constrained to a predefined 
re^onse outside a finite region in the spatial or temporal - 

5 spatial domain and in each step of the iterative. process the 
filter is constrained to a predefined tolerance in the 
wavennntber or f requency-wavenumber domain. 

6, The method of claim 1 wherein the filter is 

IQ obtained by applying alternating projection onto constraints 
defiTiing convex sets of square summable sequences « 

7 . The Kiethod^^bf claim 1 wherein the transform 

santplingyperiodicity matrix of the transform in Cairtesian 
15 coordinates is non-diagonal. 

8. nue method of claim 1, fTxcther comprising the step 
of distributing groups of receivers or single sensor seismic, 
receivers so as to obtain seismic measurements on a 

20 staggered or hexagonal grid, 

9. The method of. claim 8 wherein the step of 
transforming con^rises the use of a spatially staggered or 
hexagonal t rans format ion • 

25 

10. The it^thod of claim S "wherein the step of 
transforming the signals coinprises the use of a spatially 
staggered or hexagonal Fourier transformation. 

30 11- The method of claim 1 wherein the filter is a sero-^ 

phase finite impulse response (FIR) filter- 
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12 . Th.e method of claim 1 wherein the filter Jias at 
least two dimeasions- 

13 . a?he method of claim 1 whsreia the filter is a 3D 
5 filter. 
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ABSTRACT 

An iterative process is described for obtaining a finite 
iinpulse filter to remove noise from seismic data tlxrough an 
5 iterative proceGS with c«=nstraints on tlae filter applied in 
bofclx, the ori&inal space-tiine and the transform frequency- 
wavemairiber at each iteration step. 
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